/*
Borealis Library
A library of commonly used hydrological modeling building blocks

Copyright (C) 2008  Chris Marsh

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef RASTER_H
#define RASTER_H

#include <string>

#include <boost\shared_array.hpp>
#include <boost/noncopyable.hpp>

#include <tbb\concurrent_vector.h>

#include "VariableBase.hpp"

namespace Borealis
{
	/*
	Class: RasterPoint
		A simple X,Y,Z point for use with rasters
	*/
	template<typename zType>
	class RasterPoint
	{
	public:
		RasterPoint(unsigned int row, unsigned int col, zType z)
		{
			this->row = row;
			this->col = col;
			this->z = z;
		}
		unsigned int row;
		unsigned int col;
		zType z;
	};

	template<typename CellType>
	class RasterVisitor
	{
	public:
		virtual void operator()(CellType* cell, unsigned int row, unsigned int col) =0;
	};

	//basic contract for a raster to have
	template<typename T>
	class Raster : public VariableBase
	{
	public:
		/*
			Function: ApplyVisitor
				 Apply a vistor to each cell. Function must be independent of other cells.
			
			Parameters: 
				RasterVisitor<T> * v - Visitor object
			
			Throws:
				Never
			
			Returns:   
				void - 
		*/
		void ApplyVisitor(RasterVisitor<T>* v);

		/*
			Function: operator=
				 Assignment operator
			
			Parameters: 
				const Raster & rhs - Right hand side
			
			Throws:
				Never
			
			Returns:   
				Borealis::Raster
		*/
		Raster<T>& operator=(const Raster& rhs);


		/*
			Function: Raster
				Copy constructor.
			
			Parameters: 
				Raster & in - Raster to copy
			
			Throws:
				Never
		
			Returns:   
				 
		*/
		Raster(Raster& in);

		/*
		Function: Raster
			Default constructor

		Parameters: 
			None

		Throws:
			Never

		Returns:   

		*/
		Raster();

		/*
		Function: Raster
			Creates an empty raster of size row x col

		Parameters: 
			unsigned int row - Row count
			unsigned int col - Column count
			unsigned int bands - Raster band. Default = 0. 0 <= band < n

		Throws:
			Never

		Returns:   

		*/
		Raster(unsigned int row, unsigned int col, unsigned int bands = 0);

		/*
		Function: ~Raster
			Cleans up.

		Parameters: 

		Throws:
			Never

		Returns:   

		*/
		virtual ~Raster();

		/*
		Function: GetRow
			Returns the number of rows in the raster, 0 if empty.

		Parameters: 
			None

		Throws:
			Never

		Returns:   
			unsigned int - Row count
		*/
		inline unsigned int GetRow() const;


		/*
		Function: GetCol
			Returns the number of columns in the raster, 0 if empty.

		Parameters: 
			None

		Throws:
			Never

		Returns:   
			unsigned int - Column count.
		*/
		inline unsigned int GetCol() const;

		/*
		Function: GetRasterAt
			Returns the item at (row,col) in the raster. Used only to view the information.
			Use <SetRasterAt> to change an item.

		Parameters: 
			unsigned int row - Row 
			unsigned int col - Column
			unsigned int band - Raster band. Default = 0; 0 <= band < n

		Throws:
			std::runtime_error if a raster has not been loaded or created

		Returns:   
			T - The item  at row,col
		*/
		inline T GetRasterAt(unsigned int row, unsigned int col, unsigned int band = 0) const;


		/*
		Function: SetRasterAt
			Sets the raster to contain item 'item' at (row,col)

		Parameters: 
			unsigned int row - Row
			unsigned int col - Column	
			unsigned int band - Raster band. Default=0; 0<=band<n
			T item - Item to place

		Throws:
			std::runtime_error if a raster has not been loaded or created

		Returns:   
			void 
		*/
		inline void SetRasterAt(unsigned int row, unsigned int col,  T item, unsigned int band = 0);


		/*
		Function: CreateEmptyRaster
			Creates an empty raster. Allocates the needed memory. Not thread safe.

		Parameters: 
			unsigned int row - Row
			unsigned int col - Column
			unsigned int bands - Number of raster bands to create. Default=1. Bands must be >= 1

		Throws:
			std::bad_alloc if the memory can't be allocated

		Returns:   
			void 
		*/
		void CreateEmptyRaster(unsigned int row, unsigned int col, unsigned int bands = 1);


		/*
			Function: OpenFromFile
				Open a raster from a file
			
			Parameters: 
				std::string path - Path to file to open
			
			Throws:
				std::runtime_error if file cannot be found or error opening file
		
			Returns:   
				void 
		*/
		virtual void OpenFromFile(std::string path) = 0;

		/*
			Function: SaveToFile
				Saves the raster to a file 
			
			Parameters: 
				std::string path - Path to the file to save
			
			Throws:
				std::runtime_error if file cannot be saved
		
			Returns:   
				void  
		*/
		virtual void SaveToFile(std::string path) = 0;

		/*
			Function: GetBands
				Returns the number of bands in the raster. 
				0 = Raster is empty
				1 or more = if the raster has one or more bands
			
			Parameters: 
				None

			Throws:
				Never
		
			Returns:   
				unsigned int - Number of bands present in the raster
		*/
		virtual unsigned int GetBands() const;

		/*
		Function: GetXllCorner
		Returns x-lower-left corner extent

		Parameters: 
		None

		Throws:
		Never

		Returns:   
		double - xllcorner as defined in the raster file
		*/
		double GetXllCorner() const;

		/*
		Function: GetYllCorner
		Returns y-lower-left corner extent

		Parameters: 
		None

		Throws:
		Never

		Returns:   
		double - yllcorner as defined in the raster file
		*/
		double GetYllCorner() const;

		/*
		Function: GetCellSize
		Returns the cell size				

		Parameters: 
		None

		Throws:
		Never

		Returns:   
		unsigned int - cellsize as defined in the raster file
		*/
		unsigned int GetCellSize() const;


		/*
			Function: SetXllCorner
				 Sets the easting lower left corner (UTM).
				 Strictly speaking, this does not have to be set.
			
			Parameters: 
				double xll - UTM coordinate of easting, lower left
			
			Throws:
				Never
			
			Returns:   
				void - 
		*/
		void SetXllCorner(double xll); 


		/*
			Function: SetYllCorner
				 Sets the northing lower left corner.
				 Strictly speaking, this does not have to be set.
			
			Parameters: 
				double yll - UTM coordinate of northing, lower left
			
			Throws:
				Never
			
			Returns:   
				void - 
		*/
		void SetYllCorner(double yll);


		/*
			Function: SetCellSize
				 Set the grid cell size
			
			Parameters: 
				unsigned int cell - Cell size
			
			Throws:
				Never
			
			Returns:   
				void - 
		*/
		void SetCellSize(unsigned int cell);

	protected:
		unsigned int m_rows;
		unsigned int m_cols;
		unsigned int m_bands;

		double m_xllCorner;
		double m_yllCorner;

		unsigned int m_cellSize;

		tbb::concurrent_vector< tbb::concurrent_vector<float> > m_array;

	private:
		class ApplyVisitor
		{
			tbb::concurrent_vector< tbb::concurrent_vector<float> >* m_raster;
			RasterVisitor* m_visitor;
		public:
			ApplyVisitor(tbb::concurrent_vector< tbb::concurrent_vector<float> >* r, RasterVisitor& v)
				: m_raster(r), m_visitor(v)
			{

			}

			void operator()(const tbb::blocked_range2d<unsigned int>& r) const
			{

				tbb::concurrent_vector< tbb::concurrent_vector<float> >* raster = m_raster;
				RasterVisitor* visitor = m_visitor;

				if(!raster)
					throw std::runtime_error("raster is null");
				if(!visitor)
					throw std::runtime_error("visitor is null");

				for(unsigned int x = r.rows().begin(); x!=r.rows().end();++x)
				{
					for(unsigned int y = r.cols().begin(); y!=r.cols().end(); ++y)
					{
						visitor(raster[0][row*m_cols+col]);
					}
				}
			}
		};

	};

	template<typename T>
	void Raster<T>::ApplyVisitor( RasterVisitor<T>* v )
	{

		tbb::blocked_range2d<unsigned int> range(0,m_rows,0,m_cols);
		tbb::parallel_for(range,ApplyVisitor(m_arrary,v))
	}

	template<typename T>
	Raster<T>& Raster<T>::operator=( const Raster& rhs )
	{
// 		if(this == &rhs)
// 			return *this;



		//optimize this later

		for(unsigned int band = 0; band < m_bands; band++)
		{
			for(unsigned int row = 0; row < m_rows; row++)
			{
				for(unsigned int col = 0; col < m_cols; col++)
				{
					m_array[band].clear();
				}
			}
		}

		m_rows = rhs.GetRow();
		m_cols = rhs.GetCol();
		m_bands = rhs.GetBands();

		//allocate the # of bands
		m_array.grow_by(m_bands);

		for(unsigned int i = 0; i < m_bands; i++)
		{
			m_array[i].grow_by(m_rows*m_cols);
		}

		for(unsigned int band = 0; band < m_bands; band++)
		{
			for(unsigned int row = 0; row < m_rows; row++)
			{
				for(unsigned int col = 0; col < m_cols; col++)
				{
					m_array[band][row*m_cols+col] = rhs.GetRasterAt(row,col);
				}
			}
		}

		return *this;
	}

	template<typename T>
	unsigned int Raster<T>::GetBands() const
	{
		return m_bands;
	}

	template<typename T>
	Raster<T>::Raster( Raster<T>& in )
	{


		//optimize this later

		for(unsigned int band = 0; band < m_bands; band++)
		{
			for(unsigned int row = 0; row < m_rows; row++)
			{
				for(unsigned int col = 0; col < m_cols; col++)
				{
					m_array[band].clear();
				}
			}
		}

		m_rows = in.GetRow();
		m_cols = in.GetCol();
		m_bands = in.GetBands();
		m_xllCorner = in.GetXllCorner();
		m_yllCorner = in.GetYllCorner();
		m_cellSize = in.GetCellSize();


		//allocate the # of bands
		m_array.grow_by(m_bands);

		for(unsigned int i = 0; i < m_bands; i++)
		{
			m_array[i].grow_by(m_rows*m_cols);
		}

		for(unsigned int band = 0; band < m_bands; band++)
		{
			for(unsigned int row = 0; row < m_rows; row++)
			{
				for(unsigned int col = 0; col < m_cols; col++)
				{
					m_array[band][row*m_cols+col] = in.GetRasterAt(row,col);
				}
			}
		}
		
	}

	template<typename T>
	Raster<T>::~Raster()
	{

		
	}
	template<typename T>
	Raster<T>::Raster( unsigned int row, unsigned int col, unsigned int bands)
	{
		CreateEmptyRaster(row,col,bands);
		m_xllCorner = 0;
		m_yllCorner = 0;
		m_cellSize = 0;
	}


	template<typename T>
	void Raster<T>::CreateEmptyRaster( unsigned int row, unsigned int col, unsigned int bands )
	{
		m_rows = row;
		m_cols = col;
		m_bands = bands;

		m_array.clear();

		//allocate the # of bands
		m_array.grow_by(bands);

		for(unsigned int i = 0; i < bands; i++)
		{
			m_array[i].grow_by(row*col);
		}
	}


	template<typename T>
	void Raster<T>::SetRasterAt( unsigned int row, unsigned int col,  T item , unsigned int band)
	{

		if(m_array.empty())
		{
			throw std::runtime_error("No raster in memory");
		}

		m_array[band][row*m_cols+col] = item;
	}


	template<typename T>
	T Raster<T>::GetRasterAt( unsigned int row, unsigned int col, unsigned int band) const
	{
		if(m_array.empty())
		{
			throw std::runtime_error("No raster in memory");
		}

		return m_array[band][row*m_cols+col];
	}


	template<typename T>
	unsigned int Raster<T>::GetCol() const
	{ 
			return m_cols;
	}


	template<typename T>
	unsigned int Raster<T>::GetRow() const
	{
			return m_rows;
	}

	template<typename T>
	Raster<T>::Raster()
	{
		m_rows = 0;
		m_cols = 0;
		m_bands = 1;

	}

	template<typename T> double Raster<T>::GetXllCorner() const { return m_xllCorner;}
	template<typename T> double Raster<T>::GetYllCorner() const { return m_yllCorner;}
	template<typename T> unsigned int Raster<T>::GetCellSize() const  { return m_cellSize;}

	template<typename T> void Raster<T>::SetXllCorner(double xll) { m_xllCorner = xll;}
	template<typename T> void Raster<T>::SetYllCorner(double yll) { m_yllCorner = yll;}
	template<typename T> void Raster<T>::SetCellSize(unsigned int cell) { m_cellSize = cell;}


}
#endif